;$Id: demo.pro 4 2009-07-01 08:10:55Z dingyy.pku $
;+
; NAME: 
;       DEMO.PRO
; PURPOSE:
;       Demostration of Mie Scattering properties.
; STEPS:
;       1. Single particle properties, Dqsc vs. Dx, etc.
;       2. Lognormal distribution properties, comparision of mie_dist and mie_lognormal
;       3. Junge distribution
;-
pro demo

  ; set plot device and path
  set_plot,'z',/copy
  !P.Color=0
  !P.Background=255
  Device,set_resolution=[1500,1000]
  !P.Charsize=2
  !P.thick=2
  !P.charthick=2
  !X.Margin=[10,7]
  LOADCT,0
  TVLCT,255,0,0,1
  TVLCT,0,255,0,2
  TVLCT,0,0,255,3
  TVLCT,red,green,blue,/get
  outputPath = '~/workspace/mie_output/'
  if ~File_test(outputPath,/directory) then file_mkdir,outputPath

  ; demo of mie_single
  Dx = (dindgen(1000)+1)*0.1
  Cm = DComplex(1.42,-0.002)
  angles = dindgen(100)*180d/100d
  dqv = Cos(angles*!dpi/180d)
  mie_single,Dx,Cm,Dqv=dqv,Dqxt,Dqsc,Dqbk,Dg,Xs1,Xs2,Dph
  !P.Multi = [0,2,2,0,1]
  Plot,Dx,Dqxt,$
    xtitle='Particle size parameter',$
    ytitle='Extinction efficiency'
  Axis,YAxis=1,/Save,Ytitle='Absorption efficiency'
  OPlot,Dx,Dqxt-Dqsc,linestyle=2
  Plot,angles,Dph[*,999],/ylog,xstyle=1,$
    xtitle='Scattering angle (degree)',$
    ytitle='Phase function (X=100.0)'
  Axis,YAxis=1,YRange=[0.5,2],ylog=0,/Save,Ytitle='Phase function (X=0.1)'
  OPlot,angles,Dph[*,0],linestyle=2
  
  Plot,Dx,Dg,xtitle='Particle size parameter',ytitle='Asymmetry parameter'

  Plot,Dx,Dqbk,xtitle='Particle size parameter',ytitle='Backward scattering efficiency'

  Write_png,outputPath+'single_demo.png',TVRD(),red,green,blue

  !P.Multi = 0
  tmp_idx = 0
;  print,Xs1[*,tmp_idx]
;  print,Imaginary(Xs1[*,tmp_idx])
  Plot,angles,Imaginary(Xs1[*,tmp_idx]);,thick=1,/ylog,xstyle=1,yrange=[0,10],linestyle=2
;  OPlot,angles,Xs1[*,tmp_idx]
;  Plot,angles,Xs2[*,tmp_idx],color=3
;  OPlot,angles,Imaginary(Xs2[*,tmp_idx]),linestyle=2,color=3
;  Plot,angles,2d0 * double(Xs1[*,0]*CONJ(Xs1[*,0]) + Xs2[*,0]*CONJ(Xs2[*,0]))/(Dx[0]^2*Dqsc[0])

;  Shade_surf,Xs1,angles,Dx
  Write_png,outputPath+'single_demo_Xs1.png',TVRD(),red,green,blue

  ; Lognormal
  N = 100
  Rm = 0.1
  Sigma = 1.54
  GenLogNormal,n_r,r,N=N,d0=Rm,Sigma=Sigma
  n_logr = lognormal(N, Rm, r, Sigma) 

  n_lambda = 100
  Bext = findgen(n_lambda)
  Bsca = findgen(n_lambda)
  Kext = findgen(n_lambda)
  Ksca = findgen(n_lambda)
  Kbk  = findgen(n_lambda)
  w  = findgen(n_lambda)
  g  = findgen(n_lambda)
  Eta = findgen(n_lambda,100)
  lambda = 0.05d + findgen(n_lambda)/n_lambda * 0.7
  wavenumber = 1d/lambda
  for i = 0,n_lambda-1 do begin
    Mie_dist,n_r,r,wavenumber[i],Cm,Dqv=dqv,iKext,iKsca,iKbk,ieta
    Mie_lognormal,N,Rm,Sigma,wavenumber[i],Cm,iBext,iBsca,iw,ig
    Bext[i] = iBext
    Bsca[i] = iBsca
    w[i] = iw
    g[i] = ig
    Kext[i] = iKext
    Ksca[i] = iKsca
    Kbk[i] = iKbk
    Eta[i,*] = ieta
    print,'Step ',i,' completed.'
  endfor
  
  !P.Multi = [0,2,2,0,1]
  plot,r,n_r,/xlog,/ylog,$
    xtitle='Radius (micron)',$
    ytitle='n(r) (cm!U-3!Nmicron!U-1!N)'
  Axis,YAxis=1,/Save,Ytitle='n(log(r)) (cm!U-3!N)'
  oplot,r,n_logr,linestyle=2
;  Write_png,outputPath + 'lognormal.png',TVRD(),red,green,blue
  
  Plot,lambda,Kext,xtitle='lambda (micron)',ytitle='Extinction cross section ((10!U6!Nm)!U-1!N)'

  Axis,YAxis=1,/Save,ytitle='Scattering cross section ((10!U6!Nm)!U-1!N)'
  OPlot,lambda,Ksca,linestyle=2
;  Oplot,lambda,Bext,linestyle=2
;  Shade_surf,Alog(Eta),lambda,angles
;  Contour,Alog(Eta),lambda,angles,xtitle='Lambda (micron)',ytitle='Angle (degree)',c
;  Plot,angles,Eta[0,*]
  Plot,lambda,g,xtitle='Lambda (micron)',ytitle='Asymmetry parameter'
  Plot,lambda,w,xtitle='Lambda (micron)',ytitle='Single scatter albedo'

  Write_png,outputPath+'LogNormal_mie.png',TVRD(),red,green,blue

  ; Junge distribution
  C = 6e-5
  N_junge = 1d3
  sizeof_r = 100
  r = findgen(sizeof_r)/sizeof_r * 1 + 0.01
  dr = 1d/1000d
  nu = [4.0, 3.5, 3.0, 2.5]
  n_nu = n_elements(nu)
  Kext = fltarr(n_lambda,n_nu)
  Ksca = Kext
  Kbk = Kext
  Eta = fltarr(n_lambda,n_nu,100)
  n_r_junge = fltarr(sizeof_r,n_nu)

  for i_nu = 0, n_nu-1 do begin
    n_r = junge_r(C,nu[i_nu],r)
    n_r_junge[*,i_nu] = n_r

    tmpN = total(n_r)*dr
    print,'junge N=',tmpN
    C = C * N_junge / tmpN
    print,'Modified C=',C
    n_r = junge_r(C,nu[i_nu],r)
    print,'junge N=',total(n_r)*dr

    for i = 0, n_lambda-1 do begin
      mie_dist,n_r,r,wavenumber[i],Cm,Dqv=dqv,iKext,iKsca,iKbk,ieta
      Kext[i,i_nu] = iKext
      Ksca[i,i_nu] = iKsca
      Kbk[i,i_nu] = iKbk
      Eta[i,i_nu,*] = ieta
    endfor
  endfor

  !P.Multi = 0

  ; plot Junge distribution
;  Plot,r,n_r_junge[*,0],/xlog,/ylog,$
  Plot,r,n_r_junge[*,1],/xlog,/ylog,$
    xtitle='Radius (micron)',$
    ytitle='Junge n(r) (cm!U-3!Nmicron!U-1!N)'
  for i_nu = 2,n_nu-1 do begin
    OPlot,r,n_r_junge[*,i_nu],color=i_nu
  endfor
  Write_png,outputPath+'Junge_dist.png',TVRD(),red,green,blue
;  Plot,angles,eta[0,*,0],/ylog
;  for i=1,n_nu-1 do begin
;    OPlot,angles,eta[0,*,i

end

